************************************************************************************
*** Replication of all empirical results for Laron K. Williams, "Guns Yield Butter?  An Exploration of Defense Spending Preferences", Journal of Conflict Resolution
***
*** NOTE:
************************************************************************************

*** Set working directory
*cd ""

*** Run the do file containing necessary subprograms
do "Programs.do"

clear
clear matrix
clear mata
set seed 648
set matsize 1200
set maxvar 10000

*** Open data
use "Defense Spending Preferences.dta", clear

************************************************************************************
*** Global macros for model specifications and scenarios
************************************************************************************
global rhs_jobs "spend_health create_jobs_favor jobs_health cut_spending male age unemployed educ milex_t gdppcgrowth_tm1 us_ch_milex_tm1 alliance_us alliance_us_us mid_count_h_6"
global rhs_jobs_noint "spend_health create_jobs_favor cut_spending male age unemployed educ milex_t gdppcgrowth_tm1 us_ch_milex_tm1 alliance_us alliance_us_us mid_count_h_6"

global rhs_jobs_v1 "spend_health create_jobs_favor jobs_health cut_spending male age unemployed educ milex_pc_t gdppcgrowth_tm1 us_ch_milex_tm1 alliance_us alliance_us_us mid_count_h_6"
global rhs_jobs_v2 "spend_health create_jobs_favor jobs_health cut_spending male age unemployed educ milex_gdp_t gdppcgrowth_tm1 us_ch_milex_tm1 alliance_us alliance_us_us mid_count_h_6"

global rhs_jobs_e "spend_educ create_jobs_favor jobs_educ cut_spending male age unemployed educ milex_t gdppcgrowth_tm1 us_ch_milex_tm1 alliance_us alliance_us_us mid_count_h_6"
global rhs_jobs_p "spend_pensions create_jobs_favor jobs_pensions cut_spending male age unemployed educ milex_t gdppcgrowth_tm1 us_ch_milex_tm1 alliance_us alliance_us_us mid_count_h_6"
global rhs_jobs_u "spend_unem create_jobs_favor jobs_unem cut_spending male age unemployed educ milex_t gdppcgrowth_tm1 us_ch_milex_tm1 alliance_us alliance_us_us mid_count_h_6"

global rhs_jobs_cw "spend_health create_jobs_favor jobs_health cut_spending male age unemployed educ milex_t gdppcgrowth_tm1"
global rhs_jobs_cw_noint "spend_health create_jobs_favor cut_spending male age unemployed educ milex_t gdppcgrowth_tm1"

global rhs_jobs_cw_v1 "spend_health create_jobs_favor jobs_health cut_spending male age unemployed educ milex_pc_t gdppcgrowth_tm1"
global rhs_jobs_cw_v2 "spend_health create_jobs_favor jobs_health cut_spending male age unemployed educ milex_gdp_t gdppcgrowth_tm1"

global rhs_jobs_cw_e "spend_educ create_jobs_favor jobs_educ cut_spending male age unemployed educ milex_t gdppcgrowth_tm1"
global rhs_jobs_cw_p "spend_pensions create_jobs_favor jobs_pensions cut_spending male age unemployed educ milex_t gdppcgrowth_tm1"
global rhs_jobs_cw_u "spend_unem create_jobs_favor jobs_unem cut_spending male age unemployed educ milex_t gdppcgrowth_tm1"

*** Global macros for scenarios for substantive effects
global supportive spend_health 5 create_jobs_favor 1 jobs_health 5 cut_spending 1 male 0 age p75 unemployed 0 educ 0 milex_t p25 gdppcgrowth_tm1 p75 us_ch_milex_tm1 1.53 alliance_us 1 alliance_us_us 1.53 mid_count_h_6 0
global neither spend_health 3 create_jobs_favor 0 jobs_health 0 cut_spending 1 male 0 age p50 unemployed 0 educ 0 milex_t p50 gdppcgrowth_tm1 p50 us_ch_milex_tm1 -.52 alliance_us 1 alliance_us_us -.52 mid_count_h_6 0
global unsupportive spend_health 1 create_jobs_favor 0 jobs_health 0 cut_spending 0 male 1 age p25 unemployed 1 educ 1 milex_t p75 gdppcgrowth_tm1 p25 us_ch_milex_tm1 -5.75 alliance_us 1 alliance_us_us -5.75 mid_count_h_6 0

global supportive_us spend_health 3 create_jobs_favor 1 jobs_health 3 cut_spending 1 male 0 age p75 unemployed 0 educ 0 milex_t p25 gdppcgrowth_tm1 p75
global neither_us spend_health 3 create_jobs_favor 0 jobs_health 0 cut_spending 1 male 0 age p50 unemployed 0 educ 0 milex_t p50 gdppcgrowth_tm1 p50
global unsupportive_us spend_health 1 create_jobs_favor 0 jobs_health 0 cut_spending 0 male 1 age p25 unemployed 1 educ 1 milex_t p75 gdppcgrowth_tm1 p25

*** Include the Cold War surveys?
*keep if year >= 1988

*** Generate some survey and country fixed effects 
qui ologit spend_defense $rhs_jobs
tab id if e(sample), gen(id_)
tab ccode if e(sample), gen(cc_)

cap log close
log using "Defense Spending Preferences--Replication.smcl", replace

**********************************************************************************
*** Summary Statistics
**********************************************************************************
preserve
	qui ologit spend_defense $rhs_jobs
	keep if e(sample)
	
	qui tab ccode
	di "Countries = " r(r)
	
	qui tab id
	di "Number of elections = " r(r)
	
	sum spend_defense spend_health spend_educ spend_pensions spend_unem create_jobs_favor jobs_health cut_spending male age unemployed educ milex_t gdppcgrowth_tm1 us_ch_milex_tm1 alliance_us alliance_us_us mid_count_h_6
	
	foreach v of varlist spend_defense spend_health create_jobs_favor jobs_health cut_spending male age unemployed educ {
		tab `v'
	}
	
	collapse (count) obs = year (min) min_ts = fw_start (max) max_ts = fw_end, by(ccode studytyp)
	
	labels_ABC ccode
	list ccode studytyp min_ts max_ts obs, sepby(ccode)
restore

preserve
	tempvar more_jc more_defense more_health mn_jc mn_defense mn_health
	gen `more_jc' = cond(inlist(create_jobs_favor, .a, .b, .c), ., cond(create_jobs_favor == 1, 1, 0))
	gen `more_defense' = cond(inlist(spend_defense, .a, .b, .c), ., cond(inrange(spend_defense, 4, 5), 1, 0))
	gen `more_health' = cond(inlist(spend_health, .a, .b, .c), ., cond(inrange(spend_health, 4, 5), 1, 0))
	
	foreach v in jc defense health {
		bys id: egen `mn_`v'' = mean(100 * `more_`v'')
	}

	
	duplicates drop id, force
	foreach v in jc defense health {
		gsort -`mn_`v''
		gen rank_`v' = _n
	}
	
	labels_ABC ccode
	list ccode year rank_defense `mn_defense' rank_health `mn_health' rank_jc `mn_jc', sepby(ccode)	
restore

**********************************************************************************
*** Battery checks
**********************************************************************************

*** Correlation between spending and defense preferences
pwcorr spend_defense spend_health, sig
bys create_jobs_favor: pwcorr spend_defense spend_health, sig

*** What % have different responses across the two questions?
tab2 spend_defense_3 spend_health_3, cell

**********************************************************************************
*** Replication models
**********************************************************************************

*** Model 1: Full model with US
oprobit spend_defense $rhs_jobs cc_1 - cc_25 
qui tab id if e(sample)
di r(r)

*** Model 2: US only 
oprobit spend_defense $rhs_jobs_cw if ccode == 2
qui tab id if e(sample)
di r(r)

*** Model 3: Full model without US
oprobit spend_defense $rhs_jobs cc_2 - cc_25 if ccode != 2
qui tab id if e(sample)
di r(r)

**********************************************************************************
*** Appendix: Robustness Checks
**********************************************************************************

*** Model 4 and 5: other measures of expenditures (per capita and % of GDP)
oprobit spend_defense $rhs_jobs_v1 cc_1 - cc_25 
qui tab id if e(sample)
di r(r)

oprobit spend_defense $rhs_jobs_v2 cc_1 - cc_25 
qui tab id if e(sample)
di r(r)

*** Models 6-8: Replicate Models 1-3, but without Module I
oprobit spend_defense $rhs_jobs cc_1 - cc_25 if !inrange(year, 1985, 1986)
qui tab id if e(sample)
di r(r)

oprobit spend_defense $rhs_jobs_cw if ccode == 2 & !inrange(year, 1985, 1986)
qui tab id if e(sample)
di r(r)

oprobit spend_defense $rhs_jobs cc_2 - cc_25 if ccode != 2 & !inrange(year, 1985, 1986)
qui tab id if e(sample)
di r(r)

*** Models 9-11: Education spending preferences
oprobit spend_defense $rhs_jobs_e cc_1 - cc_25 
qui tab id if e(sample)
di r(r)

oprobit spend_defense $rhs_jobs_cw_e if ccode == 2
qui tab id if e(sample)
di r(r)

oprobit spend_defense $rhs_jobs_e cc_2 - cc_25 if ccode != 2
qui tab id if e(sample)
di r(r)

*** Models 12-14: Pensions spending preferences
oprobit spend_defense $rhs_jobs_p cc_1 - cc_25 
qui tab id if e(sample)
di r(r)

oprobit spend_defense $rhs_jobs_cw_p if ccode == 2
qui tab id if e(sample)
di r(r)

oprobit spend_defense $rhs_jobs_p cc_2 - cc_25 if ccode != 2
qui tab id if e(sample)
di r(r)

*** Models 15-17: Unemployment spending preferences
oprobit spend_defense $rhs_jobs_u cc_1 - cc_25 
qui tab id if e(sample)
di r(r)

oprobit spend_defense $rhs_jobs_cw_u if ccode == 2
qui tab id if e(sample)
di r(r)

oprobit spend_defense $rhs_jobs_u cc_2 - cc_25 if ccode != 2
qui tab id if e(sample)
di r(r)

**********************************************************************************
*** Substantive Effects
**********************************************************************************

cap drop zzz*
estsimp oprobit spend_defense $rhs_jobs cc_1 - cc_25 , genname(zzz)

*** Changes in predicted probabilites of supporting more (more + much more) defense spending
qui foreach i in supportive neither unsupportive {
	nois {
		di _newline(3) "************************************************"
		di "Scenario is `i'; Baseline probabilities"
	}
	setx $`i' 
	tempvar base4 base5 base45
	simqi, prval(4 5) genpr(`base4' `base5') level(90)
	gen `base45' = `base4' + `base5'
	qui sum `base45', meanonly
	local mn = r(mean)
	_pctile `base45', perc(2.5 5 95 97.5)
	local lo_95 = round(r(r1), .0001)
	local lo_90 = round(r(r2), .0001)
	local hi_90 = round(r(r3), .0001)
	local hi_95 = round(r(r4), .0001)
	
	nois {
		di _newline(3) "Baseline probabilities"
		di "Scenario is `i'; Pr(More) = " `mn'
		di "90% confidence interval = [" `lo_90' ", " `hi_90' "]"
		di "95% confidence interval = [" `lo_95' ", " `hi_95' "]"
	}
	
	* Create jobs = 0: Health Spending
	tempvar d4 d5 d45
	setx $`i' 	
	simqi, fd(prval(4 5) genpr(`d4' `d5')) changex(spend_health 1 5) level(90)
	gen `d45' = `d4' + `d5'
	sum `d45', meanonly
	local mn = r(mean)
	_pctile `d45', perc(2.5 5 95 97.5)
	local lo_95 = round(r(r1), .0001)
	local lo_90 = round(r(r2), .0001)
	local hi_90 = round(r(r3), .0001)
	local hi_95 = round(r(r4), .0001)
	
	nois {
		di _newline(3) "Create jobs = 0: Health Spending 1 -> 5"
		di "Scenario is `i'; Change in Pr(More) = " `mn'
		di "90% confidence interval = [" `lo_90' ", " `hi_90' "]"
		di "95% confidence interval = [" `lo_95' ", " `hi_95' "]"
	}
	
	* Create jobs = 1: Health Spending
	tempvar p4_0 p5_0 p4_1 p5_1 p45_0 p45_1 d45
	setx $`i' 
	
	setx spend_health 1 create_jobs_favor 1 jobs_health 1
	qui simqi, prval(4 5) genpr(`p4_0' `p5_0') level(90)
	gen `p45_0' = `p4_0' + `p5_0'
	
	setx spend_health 5 create_jobs_favor 1 jobs_health 5
	qui simqi, prval(4 5) genpr(`p4_1' `p5_1') level(90)
	gen `p45_1' = `p4_1' + `p5_1'
	
	gen `d45' = `p45_1' - `p45_0'
	sum `d45', meanonly
	local mn = r(mean)
	_pctile `d45', perc(2.5 5 95 97.5)
	local lo_95 = round(r(r1), .0001)
	local lo_90 = round(r(r2), .0001)
	local hi_90 = round(r(r3), .0001)
	local hi_95 = round(r(r4), .0001)
	
	nois {
		di _newline(3) "Create jobs = 1: Health Spending 1 -> 5"
		di "Scenario is `i'; Change in Pr(More) = " `mn'
		di "90% confidence interval = [" `lo_90' ", " `hi_90' "]"
		di "95% confidence interval = [" `lo_95' ", " `hi_95' "]"
	}
	
	* Cut Government Spending
	tempvar d4 d5 d45
	setx $`i' 	
	simqi, fd(prval(4 5) genpr(`d4' `d5')) changex(cut_spending 0 1) level(90)
	gen `d45' = `d4' + `d5'
	sum `d45', meanonly
	local mn = r(mean)
	_pctile `d45', perc(2.5 5 95 97.5)
	local lo_95 = round(r(r1), .0001)
	local lo_90 = round(r(r2), .0001)
	local hi_90 = round(r(r3), .0001)
	local hi_95 = round(r(r4), .0001)
	
	nois {
		di _newline(3) "Cut Government Spending 0 -> 1"
		di "Scenario is `i'; Change in Pr(More) = " `mn'
		di "90% confidence interval = [" `lo_90' ", " `hi_90' "]"
		di "95% confidence interval = [" `lo_95' ", " `hi_95' "]"
	}	
	
	* Male
	tempvar d4 d5 d45
	setx $`i' 	
	simqi, fd(prval(4 5) genpr(`d4' `d5')) changex(male 0 1) level(90)
	gen `d45' = `d4' + `d5'
	sum `d45', meanonly
	local mn = r(mean)
	_pctile `d45', perc(2.5 5 95 97.5)
	local lo_95 = round(r(r1), .0001)
	local lo_90 = round(r(r2), .0001)
	local hi_90 = round(r(r3), .0001)
	local hi_95 = round(r(r4), .0001)
	
	nois {
		di _newline(3) "Male: 0 -> 1"
		di "Scenario is `i'; Change in Pr(More) = " `mn'
		di "90% confidence interval = [" `lo_90' ", " `hi_90' "]"
		di "95% confidence interval = [" `lo_95' ", " `hi_95' "]"
	}	
	
	* Age
	tempvar d4 d5 d45
	setx $`i' 	
	simqi, fd(prval(4 5) genpr(`d4' `d5')) changex(age 46 63) level(90)
	gen `d45' = `d4' + `d5'
	sum `d45', meanonly
	local mn = r(mean)
	_pctile `d45', perc(2.5 5 95 97.5)
	local lo_95 = round(r(r1), .0001)
	local lo_90 = round(r(r2), .0001)
	local hi_90 = round(r(r3), .0001)
	local hi_95 = round(r(r4), .0001)
	
	nois {
		di _newline(3) "Age: 46 -> 63"
		di "Scenario is `i'; Change in Pr(More) = " `mn'
		di "90% confidence interval = [" `lo_90' ", " `hi_90' "]"
		di "95% confidence interval = [" `lo_95' ", " `hi_95' "]"
	}		

	* College
	tempvar d4 d5 d45
	setx $`i' 	
	simqi, fd(prval(4 5) genpr(`d4' `d5')) changex(educ 0 1) level(90)
	gen `d45' = `d4' + `d5'
	sum `d45', meanonly
	local mn = r(mean)
	_pctile `d45', perc(2.5 5 95 97.5)
	local lo_95 = round(r(r1), .0001)
	local lo_90 = round(r(r2), .0001)
	local hi_90 = round(r(r3), .0001)
	local hi_95 = round(r(r4), .0001)
	
	nois {
		di _newline(3) "College: 0 -> 1"
		di "Scenario is `i'; Change in Pr(More) = " `mn'
		di "90% confidence interval = [" `lo_90' ", " `hi_90' "]"
		di "95% confidence interval = [" `lo_95' ", " `hi_95' "]"
	}		
	
	* Military expenditures
	tempvar d4 d5 d45
	setx $`i' 	
	simqi, fd(prval(4 5) genpr(`d4' `d5')) changex(milex_t 15300000 5710000000) level(90)
	gen `d45' = `d4' + `d5'
	sum `d45', meanonly
	local mn = r(mean)
	_pctile `d45', perc(2.5 5 95 97.5)
	local lo_95 = round(r(r1), .0001)
	local lo_90 = round(r(r2), .0001)
	local hi_90 = round(r(r3), .0001)
	local hi_95 = round(r(r4), .0001)
	
	nois {
		di _newline(3) "Military expenditures: 15300000 -> 5710000000"
		di "Scenario is `i'; Change in Pr(More) = " `mn'
		di "90% confidence interval = [" `lo_90' ", " `hi_90' "]"
		di "95% confidence interval = [" `lo_95' ", " `hi_95' "]"
	}		
	
	* Real GDP Per Capita Growth
	tempvar d4 d5 d45
	setx $`i' 	
	simqi, fd(prval(4 5) genpr(`d4' `d5')) changex(gdppcgrowth_tm1 3.2 5.6) level(90)
	gen `d45' = `d4' + `d5'
	sum `d45', meanonly
	local mn = r(mean)
	_pctile `d45', perc(2.5 5 95 97.5)
	local lo_95 = round(r(r1), .0001)
	local lo_90 = round(r(r2), .0001)
	local hi_90 = round(r(r3), .0001)
	local hi_95 = round(r(r4), .0001)
	
	nois {
		di _newline(3) "Real GDP per capita growth: 3.25 -> 5.6"
		di "Scenario is `i'; Change in Pr(More) = " `mn'
		di "90% confidence interval = [" `lo_90' ", " `hi_90' "]"
		di "95% confidence interval = [" `lo_95' ", " `hi_95' "]"
	}		
	
	* Non-Alliance: US Military expenditures
	tempvar d4 d5 d45
	setx $`i' 	
	simqi, fd(prval(4 5) genpr(`d4' `d5')) changex(us_ch_milex_tm1 -1.46 2.6) level(90)
	gen `d45' = `d4' + `d5'
	sum `d45', meanonly
	local mn = r(mean)
	_pctile `d45', perc(2.5 5 95 97.5)
	local lo_95 = round(r(r1), .0001)
	local lo_90 = round(r(r2), .0001)
	local hi_90 = round(r(r3), .0001)
	local hi_95 = round(r(r4), .0001)
	
	nois {
		di _newline(3) "Non-Alliance: US military expenditures: -1.46 -> 2.6"
		di "Scenario is `i'; Change in Pr(More) = " `mn'
		di "90% confidence interval = [" `lo_90' ", " `hi_90' "]"
		di "95% confidence interval = [" `lo_95' ", " `hi_95' "]"
	}		
	
	* US Alliance = 1: US Military Expenditures
	tempvar p4_0 p5_0 p4_1 p5_1 p45_0 p45_1 d45
	setx $`i' 
	
	setx us_ch_milex_tm1 -1.46 alliance_us 1 alliance_us_us -1.46
	qui simqi, prval(4 5) genpr(`p4_0' `p5_0') level(90)
	gen `p45_0' = `p4_0' + `p5_0'
	
	setx us_ch_milex_tm1 2.6 alliance_us 1 alliance_us_us 2.6
	qui simqi, prval(4 5) genpr(`p4_1' `p5_1') level(90)
	gen `p45_1' = `p4_1' + `p5_1'
	
	gen `d45' = `p45_1' - `p45_0'
	sum `d45', meanonly
	local mn = r(mean)
	_pctile `d45', perc(2.5 5 95 97.5)
	local lo_95 = round(r(r1), .0001)
	local lo_90 = round(r(r2), .0001)
	local hi_90 = round(r(r3), .0001)
	local hi_95 = round(r(r4), .0001)
	
	nois {
		di _newline(3) "US Alliance = 1: US military expenditures: -1.46 -> 2.6"
		di "Scenario is `i'; Change in Pr(More) = " `mn'
		di "90% confidence interval = [" `lo_90' ", " `hi_90' "]"
		di "95% confidence interval = [" `lo_95' ", " `hi_95' "]"
	}	

	* Hostile MIDs
	tempvar d4 d5 d45
	setx $`i' 	
	simqi, fd(prval(4 5) genpr(`d4' `d5')) changex(mid_count_h_6 0 1) level(90)
	gen `d45' = `d4' + `d5'
	sum `d45', meanonly
	local mn = r(mean)
	_pctile `d45', perc(2.5 5 95 97.5)
	local lo_95 = round(r(r1), .0001)
	local lo_90 = round(r(r2), .0001)
	local hi_90 = round(r(r3), .0001)
	local hi_95 = round(r(r4), .0001)
	
	nois {
		di _newline(3) "Hostile MIDs (6 months prior): 0 -> 1"
		di "Scenario is `i'; Change in Pr(More) = " `mn'
		di "90% confidence interval = [" `lo_90' ", " `hi_90' "]"
		di "95% confidence interval = [" `lo_95' ", " `hi_95' "]"
	}			
}

********************************************************************************
*** Relationship between health spending and defense spending, conditional on preferences for government-financed job creation
********************************************************************************

tempname sub
postfile `sub' str20 model create outcome diff diff_lo diff_hi  /*
*/	using "Substantive Effects.dta", replace

*** Model 1: Full model with US
cap drop zzz*
qui estsimp oprobit spend_defense $rhs_jobs cc_1 - cc_25, genname(zzz)

* Create jobs = 0: Health Spending 1 -> 5
tempvar d1 d2 d3 d4 d5
setx $supportive 	
simqi, fd(prval(1 2 3 4 5) genpr(`d1' `d2' `d3' `d4' `d5')) changex(spend_health 1 5) level(90)
qui foreach c of numlist 1(1)5 {
	sum `d`c'', meanonly
	local diff`c' = r(mean)
	_pctile `d`c'', perc(5 95)
	local diff`c'_lo = round(r(r1), .0001)
	local diff`c'_hi = round(r(r2), .0001)
	
	post `sub' ("All") (0) (`c') (`diff`c'') (`diff`c'_lo') (`diff`c'_hi') 
}


* Create jobs = 1: Health Spending
tempvar p1_0 p2_0 p3_0 p4_0 p5_0 p1_1 p2_1 p3_1 p4_1 p5_1 d1 d2 d3 d4 d5
setx $supportive 
	
setx spend_health 1 create_jobs_favor 1 jobs_health 1
qui simqi, prval(1 2 3 4 5) genpr(`p1_0' `p2_0' `p3_0' `p4_0' `p5_0') level(90)
	
setx spend_health 5 create_jobs_favor 1 jobs_health 5
qui simqi, prval(1 2 3 4 5) genpr(`p1_1' `p2_1' `p3_1' `p4_1' `p5_1') level(90)
qui foreach c of numlist 1(1)5 {
	gen `d`c'' = `p`c'_1' - `p`c'_0'
	sum `d`c'', meanonly
	local diff`c' = r(mean)
	_pctile `d`c'', perc(5 95)
	local diff`c'_lo = round(r(r1), .0001)
	local diff`c'_hi = round(r(r2), .0001)
	
	post `sub' ("All") (1) (`c') (`diff`c'') (`diff`c'_lo') (`diff`c'_hi') 
}	


*** Model 2: US only 
cap drop zzz*
qui estsimp oprobit spend_defense $rhs_jobs_cw if ccode == 2, genname(zzz)

* Create jobs = 0: Health Spending 1 -> 5
tempvar d1 d2 d3 d4 d5
setx $supportive_us 	
simqi, fd(prval(1 2 3 4 5) genpr(`d1' `d2' `d3' `d4' `d5')) changex(spend_health 1 5) level(90)
qui foreach c of numlist 1(1)5 {
	sum `d`c'', meanonly
	local diff`c' = r(mean)
	_pctile `d`c'', perc(5 95)
	local diff`c'_lo = round(r(r1), .0001)
	local diff`c'_hi = round(r(r2), .0001)
	
	post `sub' ("US Only") (0) (`c') (`diff`c'') (`diff`c'_lo') (`diff`c'_hi') 
}

* Create jobs = 1: Health Spending
tempvar p1_0 p2_0 p3_0 p4_0 p5_0 p1_1 p2_1 p3_1 p4_1 p5_1 d1 d2 d3 d4 d5
setx $supportive_us
	
setx spend_health 1 create_jobs_favor 1 jobs_health 1
qui simqi, prval(1 2 3 4 5) genpr(`p1_0' `p2_0' `p3_0' `p4_0' `p5_0') level(90)
	
setx spend_health 5 create_jobs_favor 1 jobs_health 5
qui simqi, prval(1 2 3 4 5) genpr(`p1_1' `p2_1' `p3_1' `p4_1' `p5_1') level(90)
qui foreach c of numlist 1(1)5 {
	gen `d`c'' = `p`c'_1' - `p`c'_0'
	sum `d`c'', meanonly
	local diff`c' = r(mean)
	_pctile `d`c'', perc(5 95)
	local diff`c'_lo = round(r(r1), .0001)
	local diff`c'_hi = round(r(r2), .0001)

	post `sub' ("US Only") (1) (`c') (`diff`c'') (`diff`c'_lo') (`diff`c'_hi') 	
}	

*** Model 3: Full model without US
cap drop zzz*
qui estsimp oprobit spend_defense $rhs_jobs cc_2 - cc_25 if ccode != 2, genname(zzz)

* Create jobs = 0: Health Spending 1 -> 5
tempvar d1 d2 d3 d4 d5
setx $supportive	
simqi, fd(prval(1 2 3 4 5) genpr(`d1' `d2' `d3' `d4' `d5')) changex(spend_health 1 5) level(90)
qui foreach c of numlist 1(1)5 {
	sum `d`c'', meanonly
	local diff`c' = r(mean)
	_pctile `d`c'', perc(5 95)
	local diff`c'_lo = round(r(r1), .0001)
	local diff`c'_hi = round(r(r2), .0001)
	
	post `sub' ("No US") (0) (`c') (`diff`c'') (`diff`c'_lo') (`diff`c'_hi') 	
}


* Create jobs = 1: Health Spending
tempvar p1_0 p2_0 p3_0 p4_0 p5_0 p1_1 p2_1 p3_1 p4_1 p5_1 d1 d2 d3 d4 d5
setx $supportive 
	
setx spend_health 1 create_jobs_favor 1 jobs_health 1
qui simqi, prval(1 2 3 4 5) genpr(`p1_0' `p2_0' `p3_0' `p4_0' `p5_0') level(90)
	
setx spend_health 5 create_jobs_favor 1 jobs_health 5
qui simqi, prval(1 2 3 4 5) genpr(`p1_1' `p2_1' `p3_1' `p4_1' `p5_1') level(90)
qui foreach c of numlist 1(1)5 {
	gen `d`c'' = `p`c'_1' - `p`c'_0'
	sum `d`c'', meanonly
	local diff`c' = r(mean)
	_pctile `d`c'', perc(5 95)
	local diff`c'_lo = round(r(r1), .0001)
	local diff`c'_hi = round(r(r2), .0001)
	
	post `sub' ("No US") (1) (`c') (`diff`c'') (`diff`c'_lo') (`diff`c'_hi') 	
}	

postclose `sub'

preserve
	use "Substantive Effects.dta", clear
	saveold "Substantive Effects.dta", version(12) replace 
restore

*** Figure out which differences are statistically significant
preserve
	use "Substantive Effects.dta", clear
	gen sig = cond((diff_lo < 0 & diff_hi < 0) | (diff_lo > 0 & diff_hi > 0), 1, 0)
	list model create outcome sig, sepby(model create)
restore

**********************************************************************************************
*** Meta-analysis of guns versus butter tradeoff: fully-specified models (no regional dummies)
*** Estimate the change in probability of selecting more (Y=4+5), given a change in health spending from 1 to 5 (conditional on create_jobs_favor)!
*** 
*** If the variable contains non-missing values for >90% of observations, it is included in the model from the master list.
**********************************************************************************************

cap drop zzz*

levelsof id, local(C)
qui foreach i of local C {

	preserve
		keep if id == `i'
		
		qui sum ccode
		local cmn = r(mean)
		qui sum year
		local ymn = r(mean)
		nois di _newline(2) "ID = `i'; country = `cmn'; year = `ymn'"

		local vlist_`i' spend_health create_jobs_favor jobs_health
		foreach v of varlist cut_spending male age unemployed educ lower_class hh_tax_high world_worse redistribute_agree progressive_tax_agree retro_threats {
			tempvar ct total 
			egen `ct' = rownonmiss(`v')
			egen `total' = total(`ct')
			local totalnm = `total'[1]
			local perc = 100 * (`totalnm' / _N)
			
			if `perc' >= 90 {
				local vlist_`i' `vlist_`i'' `v'
			}
			else {
				local vlist`i' `vlist_`i''
			}
		}
		estsimp oprobit spend_defense `vlist_`i'', genname(zzz)
		
		* Create Jobs = 0; Health Spending = 1
		tempvar pr4_s0_c0 pr5_s0_c0 
		setx spend_health 1 create_jobs_favor 0 jobs_health 0
		simqi, prval(4 5) genpr(`pr4_s0_c0' `pr5_s0_c0')
		
		* Create Jobs = 0; Health Spending = 5
		tempvar pr4_s1_c0 pr5_s1_c0 
		setx spend_health 5 create_jobs_favor 0 jobs_health 0
		simqi, prval(4 5) genpr(`pr4_s1_c0' `pr5_s1_c0')
		
		* Create Jobs = 1; Health Spending = 1
		tempvar pr4_s0_c1 pr5_s0_c1 
		setx spend_health 1 create_jobs_favor 1 jobs_health 1
		simqi, prval(4 5) genpr(`pr4_s0_c1' `pr5_s0_c1')
		
		* Create Jobs = 1; Health Spending = 5
		tempvar pr4_s1_c1 pr5_s1_c1 
		setx spend_health 5 create_jobs_favor 1 jobs_health 5
		simqi, prval(4 5) genpr(`pr4_s1_c1' `pr5_s1_c1')
			
		* Changes in predicted probabilities
		gen d0 = (`pr4_s1_c0' + `pr5_s1_c0') - (`pr4_s0_c0' + `pr5_s0_c0')
		gen d1 = (`pr4_s1_c1' + `pr5_s1_c1') - (`pr4_s0_c1' + `pr5_s0_c1')		
		gen dd = d1 - d0

		sort id
		tempfile data`i'
		save `data`i'', replace

		drop zzz* d0 d1 dd
	restore
}

use `data1', clear
qui foreach i of local C {
	if `i' == 1 {
		di
	}
	else {
		append using `data`i''
	}
}

drop if d0 == . | id == .
keep ccode year id d0* d1* dd*

*** Generate the dataset used to produce Figure 2
preserve
	collapse (mean) d0 d1 dd id, by(ccode year)
	drop if id == .
	
	sort d0
	gen id_c0 = _n
	sort d1
	gen id_c1 = _n
	
	sort ccode year
	list ccode year id id_c0 id_c1 d0 d1 dd, sepby(ccode)
	saveold "Guns vs Butter.dta", version(12) replace 
restore

*** Summary statistics of the guns/butter effect
preserve
	*** Generate mean d0 and d1, and 90th CIs
	collapse (mean) d0 d1 dd ccode year (p5) d0_lo90 = d0 d1_lo90 = d1 dd_lo90 = dd (p95) d0_hi90 = d0 d1_hi90 = d1 dd_hi90 = dd, by(id)
	
	*** What percentage of create jobs = 0 is the gvb effect negative and statistically significant?
	tempvar d0_neg d0_negsig
	gen `d0_neg' = cond(d0 < 0, 1, 0)
	gen `d0_negsig' = 1 if (d0_lo90 < 0 & d0_hi90 < 0)
	recode `d0_negsig' (.=0)
	tab `d0_neg'
	tab `d0_negsig' if `d0_neg' == 1

	*** What percentage of create jobs = 1 is the gvb effect negative and statistically significant?
	tempvar d1_neg d1_negsig
	gen `d1_neg' = cond(d1 < 0, 1, 0)
	gen `d1_negsig' = 1 if (d1_lo90 < 0 & d1_hi90 < 0)
	recode `d1_negsig' (.=0)
	tab `d1_neg'
	tab `d1_negsig' if `d1_neg' == 1
	
	*** What percentage is the difference bigger for one than the other?  Statistically significant?
	tempvar morepos morepos_ss
	gen `morepos' = cond(d1 > d0, 1, 0)
	gen `morepos_ss' = cond(d1_lo90 > d0_hi90, 1, 0)
	tab `morepos'
	tab `morepos_ss' if `morepos' == 1
	
	*** What percentage in each quadrant (clockwise, starting with I in upper-left)?
	tempvar module
	gen `module' = 1 if d0 > 0 & d1 < 0
	replace `module' = 2 if d0 > 0 & d1 > 0
	replace `module' = 3 if d0 < 0 & d1 > 0
	replace `module' = 4 if d0 < 0 & d1 < 0
	lab def quadrant 1 "+ -" 2 "+ +" 3 "- +" 4 "- -"
	lab val `module' quadrant
	tab `module'
	
	*** Differences for the US
	list ccode year d0 d1 if ccode == 2	
restore

*** Meta-analysis of complementary attitudes?
preserve
	use "Defense Spending Preferences.dta", clear
	drop if id == .
	
	keep ccode year milex_t milex_gdp_t milex_pc_t ch_milex_tm1 us_ch_milex_tm1 ussr_ch_milex_tm1 gdppcgrowth ns_health ns_defense et_milspend* et_international* et_intpeace* et_hawk* mid_count_h - mid_count_h_36 
	duplicates drop ccode year, force
	
	sort ccode year
	tempfile data
	save `data', replace 
restore

*** Predicting the guns/butter effect (when job creation = 0)
preserve
	sort ccode year
	merge ccode year using `data'
	drop if _merge == 2
	drop _merge

	bys id: gen sim = _n

	*** Most interesting model
	global dv "d0"
	global macrorhs "milex_t gdppcgrowth ns_health ns_defense et_hawk us_ch_milex_tm1"
	
	local s = 1
	qui while `s' <= 1000 {
		reg $dv $macrorhs if sim == `s'
		if `s' == 1 {
			mat B = e(b)
			mat fit = e(r2_a), e(rmse)
		}
		else {
			mat b = e(b)
			mat B = B \ b
			mat fit = fit \ e(r2_a), e(rmse)
		}
		
		*** Displays dots to let us know that Stata is still kicking.
		if mod(`s',10) == 0 {
			nois display "." _c
			if mod(`s',1000) == 0 {
				nois display ""
			}
		}	
		local s = `s' + 1
	}
	
	clear
	matname B cons, col(7) explicit
	svmat B, names(col)
	svmat fit, names(col)
	qui foreach v of varlist _all {
	
		if "`v'" == "milex_t" {
			sum `v', meanonly
			local mn = r(mean)
		
			_pctile `v', perc(2.5 5 95 97.5)
			local lo95 = r(r1)
			local lo90 = r(r2)
			local hi90 = r(r3)
			local hi95 = r(r4)
		
			nois di _newline(1) "`v': mean = " `mn' "   90% CI = [" `lo90' ", " `hi90' "] 95% CI = [" `lo95' ", " `hi95' "]"
		
		}
		else {
			sum `v', meanonly
			local mn = round(r(mean), .001)
		
			_pctile `v', perc(2.5 5 95 97.5)
			local lo95 = round(r(r1), .001)
			local lo90 = round(r(r2), .001)
			local hi90 = round(r(r3), .001)
			local hi95 = round(r(r4), .001)
		
			nois di _newline(1) "`v': mean = " `mn' "   90% CI = [" `lo90' ", " `hi90' "] 95% CI = [" `lo95' ", " `hi95' "]"
		}
	}
restore

*** Predicting the guns/butter effect (when job creation = 1)
preserve
	sort ccode year
	merge ccode year using `data'
	drop if _merge == 2
	drop _merge

	bys id: gen sim = _n

	*** Most interesting model
	global dv "d1"
	global macrorhs "milex_t gdppcgrowth ns_health ns_defense et_hawk us_ch_milex_tm1"
	
	local s = 1
	qui while `s' <= 1000 {
		reg $dv $macrorhs if sim == `s'
		if `s' == 1 {
			mat B = e(b)
			mat fit = e(r2_a), e(rmse)
		}
		else {
			mat b = e(b)
			mat B = B \ b
			mat fit = fit \ e(r2_a), e(rmse)
		}
		
		*** Displays dots to let us know that Stata is still kicking.
		if mod(`s',10) == 0 {
			nois display "." _c
			if mod(`s',1000) == 0 {
				nois display ""
			}
		}	
		local s = `s' + 1
	}
	
	clear
	matname B cons, col(7) explicit
	svmat B, names(col)
	svmat fit, names(col)
	qui foreach v of varlist _all {
	
		if "`v'" == "milex_t" {
			sum `v', meanonly
			local mn = r(mean)
		
			_pctile `v', perc(2.5 5 95 97.5)
			local lo95 = r(r1)
			local lo90 = r(r2)
			local hi90 = r(r3)
			local hi95 = r(r4)
		
			nois di _newline(1) "`v': mean = " `mn' "   90% CI = [" `lo90' ", " `hi90' "] 95% CI = [" `lo95' ", " `hi95' "]"
		
		}
		else {
			sum `v', meanonly
			local mn = round(r(mean), .001)
		
			_pctile `v', perc(2.5 5 95 97.5)
			local lo95 = round(r(r1), .001)
			local lo90 = round(r(r2), .001)
			local hi90 = round(r(r3), .001)
			local hi95 = round(r(r4), .001)
		
			nois di _newline(1) "`v': mean = " `mn' "   90% CI = [" `lo90' ", " `hi90' "] 95% CI = [" `lo95' ", " `hi95' "]"
		}
	}
restore


log close
